{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import copy\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import torch\n",
    "\n",
    "import cooper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Maximum entropy\n",
    "\n",
    "This example is inspired by [this StackExchange question](https://datascience.stackexchange.com/questions/107366/how-do-you-solve-strictly-constrained-optimization-problems-with-pytorch).\n",
    "\n",
    ">I am trying to solve the following problem using pytorch: given a six sided die whose\n",
    ">average roll is known to be 4.5, what is the maximum entropy distribution for the faces?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MaximumEntropy(cooper.ConstrainedMinimizationProblem):\n",
    "    def __init__(self, mean_constraint):\n",
    "        self.mean_constraint = mean_constraint\n",
    "        super().__init__(is_constrained=True)\n",
    "\n",
    "    def closure(self, probs):\n",
    "        # Verify domain of definition of the functions\n",
    "        assert torch.all(probs >= 0)\n",
    "\n",
    "        # Negative signed removed since we want to *maximize* the entropy\n",
    "        entropy = torch.sum(probs * torch.log(probs))\n",
    "\n",
    "        # Entries of p >= 0 (equiv. -p <= 0)\n",
    "        ineq_defect = -probs\n",
    "\n",
    "        # Equality constraints for proper normalization and mean constraint\n",
    "        mean = torch.sum(torch.tensor(range(1, len(probs) + 1)) * probs)\n",
    "        eq_defect = torch.stack([torch.sum(probs) - 1, mean - self.mean_constraint])\n",
    "\n",
    "        return cooper.CMPState(\n",
    "            loss=entropy, eq_defect=eq_defect, ineq_defect=ineq_defect\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmp = MaximumEntropy(mean_constraint=4.5)\n",
    "formulation = cooper.LagrangianFormulation(cmp)\n",
    "\n",
    "probs = torch.nn.Parameter(torch.rand(6))\n",
    "\n",
    "primal_optimizer = cooper.optim.ExtraSGD([probs], lr=3e-2, momentum=0.7)\n",
    "dual_optimizer = cooper.optim.partial_optimizer(\n",
    "    cooper.optim.ExtraSGD, lr=9e-3, momentum=0.7\n",
    "    )\n",
    "\n",
    "coop = cooper.ConstrainedOptimizer(\n",
    "    formulation=formulation,\n",
    "    primal_optimizer=primal_optimizer,\n",
    "    dual_optimizer=dual_optimizer,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CMPState(loss=tensor(-1.6136, grad_fn=<SumBackward0>), ineq_defect=tensor([-0.0544, -0.0788, -0.1142, -0.1654, -0.2398, -0.3475],\n",
      "       grad_fn=<NegBackward0>), eq_defect=tensor([3.6955e-06, 0.0000e+00], grad_fn=<StackBackward0>), proxy_ineq_defect=None, proxy_eq_defect=None, misc=None)\n"
     ]
    }
   ],
   "source": [
    "state_history = cooper.StateLogger(save_metrics=[\"loss\", \"eq_defect\", \"eq_multipliers\"])\n",
    "\n",
    "for iter_num in range(5000):\n",
    "\n",
    "    coop.zero_grad()\n",
    "    lagrangian = formulation.composite_objective(cmp.closure, probs)\n",
    "    formulation.custom_backward(lagrangian)\n",
    "    coop.step(cmp.closure, probs)\n",
    "\n",
    "    # Store optimization metrics at each step\n",
    "    partial_dict = {\"params\": copy.deepcopy(probs.data)}\n",
    "    state_history.store_metrics(formulation, iter_num, partial_dict)\n",
    "\n",
    "print(cmp.state)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(1.6136)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# True solution for mean=4.5\n",
    "a = torch.tensor([0.05435, 0.07877, 0.1142, 0.1654, 0.2398, 0.3475])\n",
    "-torch.sum(a * torch.log(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3IAAADWCAYAAACUsrCEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABoiElEQVR4nO3dd3hb5dn48e8t2ZblveOdvTeEAGGXvTeFlrYUKIW3tH27aGl/b0t3ebv7Qgu0tLRAaYEyyw5lQ4AQQvZOHI94b8uyLen5/XGk2HHkLUu2dH+uS5eko6NznhMlJ+c+z/PctxhjUEoppZRSSik1edgi3QCllFJKKaWUUiOjgZxSSimllFJKTTIayCmllFJKKaXUJKOBnFJKKaWUUkpNMhrIKaWUUkoppdQko4GcUkoppZRSSk0yGsipkBIRIyKzBvl8s4icPMxt7ROR0/yvvy0ifwpNK5VS6lAicpOI1IhIu4hkR7o9SqnoJCK3icgDg3w+7OukEe53XLarIksDOXWQP3DqFpGcfss/9Ado00a4vftE5Ed9lxljFhpjXh1p24wxPzHGXD/S7ymlYoP//NUpIm0i0iwib4vIjSIy5P9zIhIP/Ao4wxiTYoxpGGUbpvnPlXGj+b5SKjqIyDUislFEXCJSLSJ/EJGM4Xx3tNdJ/fYfsusvNbFpIKf62wtcFXgjIouBpMg1Z2zEon/PlYoN5xtjUoGpwM+AbwL3DuN7U4BEYPM4tk0pFQNE5GvA7cA3gHTgGKxz0ksikhDJtqnooxe4qr/7gU/3ef8Z4G+BNyLyqohc3+f9NSLyZv+NiMgNwCeBW/xDlZ72L+87XPI2EXlURP7pv4u+TkSWBmtU/6EIInKM/457s4h81He4gL+NPxaRtwAXMMPfzj3+/ewVkU+O6k9HKTXhGWNajDFPAR8HPiMii0TEISK/EJH9/iGUd4mIU0TmANv9X20Wkf8AiMg8EXlJRBpFZLuIXBHYvv97vxSRMhFpEZE3RcQJvN5nO+0icqyIzBKR1/zr1YvIP8P5Z6GUCh8RSQO+D3zRGPO8MabHGLMPuAKYBlztXzVxoGufftdJNhH5lojsFpEGEXlYRLL6rHt8n2uhcv+1zqDXXyJS6B+90Hc7y/3np3j/+2tFZKuINInICyIydTz/3NToaSCn+lsDpInIfBGxA1cCA47lHogx5h7gQeB//UOVzh9g1QuBR4As4O/AE4ETyUBEpAh4BviR/3tfB/4lIrl9VvsUcAOQCtQBvwPO9t+tXwWsH+kxKaUmF2PMe0AFcAJWD90cYBkwCygCvmuM2QEs9H8lwxjzMRFJBl7COiflYZ0Hfy8iC/zr/QI4EutckgXcAviAE/tsJ8UY8w7wQ+BFIBMoBv5v3A5YKRVpq7B69x/ru9AY0w48C5zuXzTca58vAhcBJwGFQBNwJ4A/uHoO65ySi3VuWz/U9Zcxpgp4B7i0z+JPAI8aY3pE5ELg28Al/u2+ATw0wj8HFSYayKlgAr1ypwNbgcpx3NcHxphHjTE9WHNUErGGIQzmauBZY8yzxhifMeYlYC1wTp917jPGbDbGeAAP1kXWIhFxGmMOGGN0CJVSsaEK62LpBuArxphGY0wb8BOsAC2Y84B9xpi/GGM8xpgPgX8Bl/uHal8LfNkYU2mM8Rpj3jbGdA2wrR6sYVWFxhi3MeawEQxKqaiRA9T7rz36O+D/HIZ/7XMj8B1jTIX/HHMbcJl/Hu4ngNXGmIf8PX8Nxpj1w2zn3/FPoxERwToX/r3PPn9qjNnqP46fAMu0V25i0kBOBXM/1gniGvoMqxwn5YEXxhgf1t3zwiG+MxXrgqo58ACOBwoG2G4H1hCrG4EDIvKMiMwLUfuVUhNbERCHNdf3gz7njOex7jYHMxU4ut855pNAPtaFWCKwe5j7vwUQ4D2xssZdO+ojUUpNdPVAjgRPeFTg/xyGf+0zFXi8z3loK+DFmtdbwvDPQ/39CzhWRAqwRhL4sHreAvv8bZ99NmKdw4pGuS81jjSQU4cxxpRhJT05h37DA4AODk1+kj/Ypoaxu5LAC/+d7mKsO+iDKQfuN8Zk9HkkG2N+NtC+jTEvGGNOxzqRbgP+OIy2KaUmMRE5Cuvi4wmgE1jY55yRboxJGeCr5cBr/c4xKcaYm7AuxNzAzCDfO+ycZ4ypNsZ8zhhTCHwea4jmgCValFKT2jtAF9awxINEJAU4G3jZv2i41z7lWNNC+p6LEo0xlf7Pgp2HYIjrL2NME9aQ749j3bj/hzEm8J1y4PP99uk0xrw92DZVZGggpwZyHfAxf29WX+uBS0QkyX8xct0g26gBZgyxnyNF5BL/3av/xjoBrhniOw8A54vImSJiF5FEETlZRIqDrSwiU0TkQv+8ly6gHevuk1IqColImoicB/wDeMAY8xHWzZtfi0ief50iETlzgE38G5gjIp8SkXj/4ygRme+/e/5n4Ff+pAF2f1ITB9Z8XB99znsicnmfc1MT1gWWnn+UikLGmBasZCf/JyJn+c8d04CHsXrd7vevOtxrn7uAHweGNYpIrn8OG1jz4E4TkStEJE5EskVkmf+z4Vx//R1rGs1l9A6rDOzzVhFZ6N9nuohcPrw/ARVuGsipoIwxu40xa4N89GugG+sk8VesE8lA7gUW+LvnnxhgnSex7gg1YSUoucQ/ZnywtpVjTRT+NtaFUzlWmt+B/j7bgK9i3e1qxJo0fNNg+1BKTUpPi0gb1jnhO1hzTz7r/+ybwC5gjYi0AquBucE24p9DdwbWvJEqoBornbjDv8rXgY3A+1jnlNsBmzHGBfwYeMt/3jsGOAp4V0Tagaew5tbtCelRK6UmDGPM/2Jdn/wCaAXexTonndpnLu1wr31+i3XeeNF/blsDHO3fz36skVNfwzoPrQcC2S+Hc/31FDAbqPbf7Aq0/3Gsc9o//OfKTVi9iWoCkt6eVKXCS0RuA2YZY64eal2llFJKqWgnIvuBq40xrw+5sop52iOnlFJKKaVUhPnLKOUC+yLcFDVJaCCnlFJKKaVUBPmTM+0E/s8/bFKpIenQSqWUUkoppZSaZLRHTimllFJKKaUmGQ3klFJKKaWUUmqSCVZ5fkLIyckx06ZNi3QzlFIh9sEHH9QbY3Ij3Y6x0POTUtEnXOcmf02u24D5wMoBSv0gIhnAn4BFWPUHrzXGvDPYtvXcpFT0GezcNGEDuWnTprF2bdBzm1JqEhORski3Yaz0/KRU9AnjuWkTcAlw9xDr/RZ43hhzmYgkAElDbVjPTUpFn8HOTRM2kFNKKaWUijbGmK0AIjLgOiKSDpwIXOP/TjfQHYbmKaUmEZ0jp5RSSik1sUwH6oC/iMiHIvInEUmOdKOUUhOLBnJKKaWUUiEkIqtFZFOQx4XD3EQccATwB2PMcqAD+NYA+7pBRNaKyNq6uroQHYFSajLQoZVKKaWUUiFkjDltjJuoACqMMe/63z/KAIGcMeYe4B6AFStWaHFgpWKIBnJKRSF3j5fGjm5c3R46u324uj24ery4u714jcHrM/iMwecDrzH4fAafsV5jgl8HDHZ1cO7iArJTHONzMONMREqAvwFTsA7zHmPMbyPbqhHydEFnM7ibrees6ZCSF+FGKaVGyxhTLSLlIjLXGLMdOBXYEul2KRXNvD5Dm7uHZlcPLZ09dHl8dHm8dHt8dHl8dHt89Hh9vddDBgwGY6yLB9PvPcYMeu109dFTsdkGnis7HBrIKTUJGWOobO5kU2Ure+s72FPXzr6GDurauqhv76a9yxPW9hw5NXPSBnKAB/iaMWadiKQCH4jIS8aYiXPR5PNBwy7obOoN1tzN/vct0OM6dP30Ejjx6xFoqFJqKCJyMfB/QC7wjIisN8acKSKFwJ+MMef4V/0i8KA/Y+Ue4LORabFS0cPd42VbdRubKlvYU9dBeZOL8kYXB1rctLp7BrqXPS4+efTUMW9DAzmlJomKJhf/2VbLmzvr+bC8mbq2roOf5aY6mJ6dzOLiDHJSEshJcZCVnECyI46keDvOBOuRGGcn3i6ICHabYBdBBOu1zXptGyyT2gDL053xIT7a8DHGHAAO+F+3ichWoIiJcvfbGFh7L9Rs6l2WkALODEjKhuxZ1uvEDHBmQv1O2PEcNJdDRkmEGq2UGogx5nHg8SDLq4Bz+rxfD6wIX8uUij4er491+5t5bUctb+ysZ0tVKx6fFa054+2UZDkpyUxixbRMspISyEhKICMpnnRnPInxdhLibDjibCTE2Uiw24i396YXEbGyz0rgNeJ/Bvq/D2KMnXGABnJKTWh1bV08/mEFj62rZFt1GwAlWU5OmJXD8tIMlhRnMCM3mdTEyRtITSQiMg1YDrw7xKrhs+0ZK4ibey4UHQmJaWAf5PdOLYBdq2H/OxrIKaWUikkbK1r417oKnv6oioaObuw24YjSDD5/0gwWF6WzqCidogznoGVAJgMN5JSagD4qb+au13bz0pYaPD7DEaUZ/L9z5/OxeXnMyE2JdPOikoikAP8C/tsY0xrk8xuAGwBKS0vD06iKD2DXS1C6Cmafbt3yG0pCEhQug8oPYMGFEDdph7wqpZRSw+bzGVZvreGe1/ewtqyJBLuN0xdM4fylBayalUNaFN701kBOqQlk3f4mfvXiDt7cVU9aYhzXHT+dy1eUMCtPg7fxJCLxWEHcg8aYx4KtE/bMcE1l8NFDkDUTFl06vCAuYOoqqHgfqj6E0mPGr41KKaXUBPDajjp++uxWtlW3UZzp5HvnL+CS5cWkJ0Vf8NaXBnJKTQBVzZ3c/vw2nlxfRU6Kg1vPnscnj5lKikP/iY43scZV3AtsNcb8KtLtAaxkJu//CRypsOJasI/w70HmdEjJh7K3NJBTSikVtfbVd/Ddpzbz+o46pmYn8dsrl3Hu4gLi7LFRKluvEpWKIJ/PcP+aMn723Da8xnDzKbO46eSZJGsAF07HAZ8CNorIev+ybxtjno1Ia7w9VhDn6YLj/wsco+iNFYGpx8Lmx6GlAtKLQ99OpZRSKkJ8PsNf39nH7c9vI95m4/+dO59PHTsVR5w90k0LK71aVCpCKppc3PLoBt7e3cCJc3L58UWLKMlKinSzYo4x5k0GTioVXsbA+r9bwddR10Fawei3VXwUbH0a9q+BxZeFro1KKaVUBNW1dfGlhz7knT0NnDI3l59esoT89MRINysiNJBTKgJe2VbLl//xIT4DP7tkMR8/qmTSZ05SIbBrNVStg3nnQ/7isW0rIRkKlllz5eafr0lPlFJKTXoflDXyXw+uo6Wzh9svXcwVK2L7+kkDOaXCyOsz/Pblnfzu5Z0sKEjjrquPpDRbe+EUUL0Rtv0bilbArFNDs82pq6ByLVSth9KjQ7NNpZRSKgIeXlvOtx/bSFGmk/s+u5L5BWmRblLEaSCnVJi4e7x8+R8f8sLmGi49opgfX7yIxPjYGsutBtBSCevuh4xSWHrlyDJUDiZrBqRMgf1vayCnlFJqUjLGcMd/dvHLl3Zwwuwc7vjEEaQ7ozsb5XBpIKdUGDR1dHP939aybn8T/3PeAq49blpMDwVQfXS1WclN4hPhqOsHL/Y9UiJQeixseQJaqyCtMHTbVkoppcaZz2f47lObeGDNfi5ZXsTPLl1CQlxsZKQcjjH/SYhIiYi8IiJbRGSziHw5yDoiIr8TkV0iskFEjhjrfpWaLKqaO7nsrrfZWNnCnZ84guuOn65BnLJ4PbD2L1Ywd9T1kJge+n2UrARbHJS9HfptK6WUUuPE5zN8+/GNPLBmP58/aQa/vGKpBnH9hOJPwwN8zRizADgG+IKILOi3ztnAbP/jBuAPIdivUhNeZXMnV96zhtrWLh647mjOWTyGLIQquhgDGx+Bxt2w7BPWsMrxkJAMBUuh8gPwdI/PPpRSSqkQ8vkM/+/JTfzj/XJuPmUW3zprnt4ED2LMgZwx5oAxZp3/dRuwFSjqt9qFwN+MZQ2QISJ6RauiWkWTiyvveYcmVzf3X380K6dnRbpJaiLZ+xqUr4HZZ0DROA9SKF0FPS448NH47kcppZQaI2MMP/j3Fv7+7n5uOnkmXztjjgZxAwhp/6SITAOWA+/2+6gIKO/zvoLDgz2losaBlk6u+uMaml09PHDd0SwryYh0k9REUrsNNj8B+Utg7jnjv7/smZCcB2Vvjf++lFJKqTG4+/U93Pf2Pq47fjq3nDlXg7hBhCyQE5EU4F/AfxtjWke5jRtEZK2IrK2rqwtV05QKq2ZXN5++9z2aOqwgbqkGcaqv9lr44D4r8cjyq0OXoXIwIjD1WGjaC60Hxn9/Siml1Cg8ub6Snz23jfOXFvKdc+ZrEDeEkARyIhKPFcQ9aIx5LMgqlUBJn/fF/mWHMMbcY4xZYYxZkZubG4qmKRVWnd1erv/rWsoaXNzz6SM1iFOH6nbBe38Em91KbhLOIt3F/qQn+98J3z6VUkqpYXp7dz1ff+QjjpmRxS8uX4LNpkHcUEKRtVKAe4GtxphfDbDaU8Cn/dkrjwFajDF6W1hFFY/XxxcfWscH+5v4zZXLWDUzJ9JNUhOJz2f1xLka4KjrICnMcyYdKdZQzor3wdsT3n0rpZRSgyhvdPFfD65jWnYyd39qBY44rbM7HKGoI3cc8Clgo4is9y/7NlAKYIy5C3gWOAfYBbiAz4Zgv0pNKN9/egurt9bywwsXanZKdbgtT0D9dlh6lVWoOxJKVkLVOmjcA7lzI9MGpZRSqg9Xt4fP/W0tPp/hT59ZocW+R2DMgZwx5k1g0L5PY4wBvjDWfSk1Ud3/zj7uX1PG50+cwaeOnRbp5qiJpuwdK0vljJOh9JjItSO92HpurdJATqkIEZHLgduA+cBKY8zaIOvMBf7ZZ9EM4LvGmN+Eo41KhYsxhlse3cD2mjb+cs1RTM1OjnSTJpVQ9MgpFdPe3lXPbU9v4WPz8rjlrHmRbo6aaBp2W/XicufD/Asj2xZHKiSkQFt1ZNuhVGzbBFwC3D3QCsaY7cAyABGxY+UVeDwcjVMqnP74xh7+veEAt5w1l5Pn5kW6OZOOBnJKjcG++g5uenAdM3KS+e2Vy7DrxNzY4O2BlsPyNR3O44a1f4akbDji02ALacWX0UnNh3YN5JSKFGPMVmAk2fhOBXYbY8rGrVFKRcAHZU3c/vx2zl6Uz00nzYx0cyYlDeSUGqX2Lg/X/20tInDvZ44iNVHHdMeM9mp4/X+Ht258Eqy8ARKSxrdNw5VaYCU8MSY8pQ+UUmN1JfBQpBuhVCi1dPbwpYc+pDAjkdsvW6JlBkZJAzmlRsEYw62PbWRPXTsPXHc0pdkT5CJdhUdSDqy4dnjrppeEP0PlYFILrJ7CzqaJ1S6looiIrAbyg3z0HWPMkyPYTgJwAXDrIOvcANwAUFpaOsKWRp+KJhf7G1wcOzNbg4MJyhjDt/61gZpWN4/ceCxpeiN81DSQU2oUHlhTxtMfVfGNM+eyapaWGYg58U4oWBrpVoxOqv/asq1aAzmlxokx5rQQbepsYJ0xpmaQfd0D3AOwYsUKE6L9Tlo/fW4bz2w4wMLCNP77tDmcNj9PA7oJ5u/v7ee5TdV86+x5LC/NjHRzJrUJMGFDqcllQ0UzP/z3Vk6Zm6tjutXkczCQ01KeSk0CV6HDKkekrKGDadlJtHdZKe2v/+thSUFVBO2oaeMHT2/hhNk53HBChErxRBEN5JQagRZXD//14DpyUx386opl2DS5iZpsEpLBkaaZK5WKEBG5WEQqgGOBZ0TkBf/yQhF5ts96ycDpwGORaenkVNHUyfGzc3j5qyfxiaNLeXlbLa5uT6SbpYAer4+v/HM9KY44vYYKEQ3klBomYwxfe2Q9Na1u7vzkEWQmJ0S6SUqNTmqB9sgpFSHGmMeNMcXGGIcxZoox5kz/8ipjzDl91uswxmQbY1oi19rJpc3dQ7Orh+LMJOLsNo6dkQ3A/kZXhFumAO74zy42V7Xyk0sWk5vqiHRzooIGckoN0wNryli9tZZvnzOfZSUZkW6OUqOXmg/tNVbmSqWUihKVzZ0AFGc6AZjqT0RW1qCBXKRtrGjhjld2ccnyIs5cGCwPkBoNDeSUGoZ99R385NltnDQnl2tWTYt0c5Qam9QC8HaDqzHSLVFKqZCpaAwEclYANzUrGbDmzanIcfd4+doj68lJSeB75y+MdHOiigZySg3B6zN87ZGPiLcLt1+qtU5UFDiY8KQqsu1QSqkQqmiyet4CPXLpSfGkO+O1Ry7Cfr16Bztq2rn90iWkJ2mpgVDSQE6pIdzz+h4+KGvihxctIj89MdLNUWrsUgusZ014opSKIhVNnSTG28juM4d9anaSzpGLoA/KGrnn9T1ctbKUk+fmRbo5UUcDOaUGsa26lV+/tIOzF+VzwdLCSDdHqdCIT4TEDA3klFJRpaKpk+LMpENGzpRmJWmPXIR0ebzc8ugGCtOdfOfc+ZFuTlTSQE6pAXR7fHzlnx+R5ozjRxct0iGVKrpo5kqlVJSpaHYdHFYZMC07mcrmTnq8vgi1Knb94dXd7K7r4McXLyLFERfp5kQlDeSUGsD//WcnWw+08tNLlpCdomlyVZRJK7AyV/r04kYpFR0qmzoPC+RKs5Pw+gxV/oyWKjx21bbx+1d2c+GyQh1SOY40kFMqiPXlzfz+1d1cdmQxpy+YEunmKBV6Kfng84CrPtItUUqpMWvv8tDkryHX19QsLUEQbj6f4Vv/2kiSw87/nLcg0s2JahrIKdWPu8fLVx9eT35aIt89X09AKkodzFypwyuVUpNfZdOhNeQCpmZrCYJwe+j9/awta+I758wnR0c0jSsN5JTq5/bnt7GnroOfX7aEtERNk6ui1MFAThOeKKUmv0DpgaKMQwO5vFQHjjib9siFSU2rm589u41VM7O57MjiSDcn6mkgp1Qfb++u5y9v7eOaVdNYNSsn0s1RYSIifxaRWhHZFOm2hE2cA5KytUdOKRUVKpoOLQYeYLMJU7OTKNMSBGFx21Ob6fb6+MnFizVJXBiEJJAb6iJIRE4WkRYRWe9/fDcU+1UqlNrcPXzjkQ3MyEnmm2fNi3RzVHjdB5wV6UaEXWqB9sgppaJCRZMLR5yNnJSEwz4rzUpmv/bIjbsXN1fz3KZqvnzabKblJEe6OTEhVD1y9zH0RdAbxphl/scPQrRfpULmh//ewoGWTn5xxVKcCfZIN0eFkTHmdaAx0u0Iu9R8K3Ol1xPpliil1JhU+DNWBusFChQFN8ZEoGWxoc3dw3ef3My8/FQ+d8KMSDcnZoQkkIvZiyAVNVZvqeHhtRXcdPJMjijNjHRz1AQkIjeIyFoRWVtXVxfp5oRGagEYH3REyfEopWJWoBh4MFOzk+js8VLX1hXmVsWOn7+wnZo2Nz+7dAnxdp25FS7hrM53rIh8BFQBXzfGbA7jvpUaUGNHN996bCPzC9L48qlzIt0cNUEZY+4B7gFYsWJF2G7rdnm87Kt3UdPqHvZ3puckU5IV/ILmEH0zV6YVjLKFSikVeRVNLpYUpwf9rDRQgqDRRV5aYjibFRM+KGvi/jVlXLNqGstKMiLdnJgSrkBuHTDVGNMuIucATwCz+68kIjcANwCUlpaGqWkqlhlj+J8nNtHS2c39160kIU7vIqnI8vkMlc2d7KprZ1dNO/sbXXh8I48bV07P5KyFBYMPE06ZAog1vFIppSapgWrIBfSWIHBx1LSscDYt6vV4fXz7sY0UpCXytTPmRro5MScsgZwxprXP62dF5PcikmOMqe+3XkTueKvY9dRHVTyz8QDfOHMu8wvSIt0cFaMaO7rZVdvOzto29tR14Or2AlCQnsiqmdnMykuhODMJ2zDuM3h9hle31/Hmrnq2HmjjvCUFLC5KD549zB4PyTnQWhXiI1JKqfAZqIZcQFGGE7tNtJbcOLj3zb1sr2njj59eQYojnAP9FIQpkBORfKDGGGNEZCXW3LyGcOxbqYHUtLr57pObWV6awedP1Im5sUxEHgJOBnJEpAL4njHm3vHaX2e3l9117eyqtR4NHd0ApDnjmF+Qxqy8FGbmJpM6yjqG5ywuYGlJBk98WMlD75WzbkoTFy4rIjP58GxupOZr5kql1KQWqCE3UCCXEGejMCNRa8mFWHmji9+s3sHpC6Zw+oIpkW5OTApJIBfsIgiIBzDG3AVcBtwkIh6gE7jSaOogFUHGGG55dANdHi+/umIZcToxN6YZY64az+17vD72N7r8vW7tVDZ3Ygw44mzMyE0+2OuWm+oIWd2dogwnN500k3f2NPDSlhp+s3oHpy2YwnEzc7DZ+uwjtQBqNoO3x+qhCyefD5r3QVIOJGqPuIoNInI5cBswH1hpjFk7wHpfAa4HDLAR+KwxZviTZWPIQDXk+irKcHKgpTNcTYp6xhi+99RmbCJ8/4KFkW5OzApJIDfURZAx5g7gjlDsS6lQeOi9cl7bUcf3L1jIdK11okaopbOHZzYMXUjbYKhv62Jfg4sujw8RKMlM4pS5eczKS6E0Kwm7bfwKptpswnGzclhYmMZTH1Xx7MZq1u9v5prjpvX29qXmW5kr22shvWjc2nKY+l2w7j7oaoP8JXDUdeHbt1KRtQm4BLh7oBVEpAj4ErDAGNMpIg8DV2KVe1L9DFZDLiAvNZGPKprD16go9/ymav6zrZb/d+58CjOC94Sq8aeDWVXM2d/g4kfPbOG4Wdl86pipkW6OmoQ6ujy8v294FVdSE+NYXprBzNwUZuamRKRGYUZSAp86Ziqbq1p5eG05j39YyaeOmWr1/qX6s1W2VYc3kNv3BhgDmdOhaa/1OkS9kUpNZMaYrcBwet/jAKeI9ABJWFm/VRCD1ZALyEt1UNvahTEmZCMfYlWbu4fbnt7M/II0rlk1LdLNiWkayKmY4vUZvv7IR9hF+PllSw8dYqbUMBVmOLltkg0lEREWFaXT0tnDvzcc4P19TaycngXJeSA2qwRBODXvh+xZ1mPTo+BuBqfWcFQKwBhTKSK/APZjTUl50RjzYrB1NeP34DXkAvLSHHT2eGnv8ox6/rGy/PLFHdS2dXHX1Ufq1JQI0z99FVP+/OZe3tvXyPcuWKhDAVRMWjUzm5m5yTy78QAN7V1gj7OCuXAGcl3t0NkIGaXWA6CpLHz7V2qcichqEdkU5HHhML+fCVwITAcKgWQRuTrYusaYe4wxK4wxK3Jzc0N3EJNIRZNrwEQnAXmpVv24Wi0KPiYbK1r42zv7+OTRpSwv1ZtvkaaBnIoZO2ra+PmL2zl9wRQuPSKMQ8iUmkBEhMuPLEEEHvmgAp/PQOqU8GaubCm3njNKIa0IbHFWD51SUcIYc5oxZlGQx5PD3MRpwF5jTJ0xpgd4DFg1fi2evIaqIReQl+oAoLZVA7nR8nh9fPvxjWSnOPjGmfMi3RyFBnIqRvR4fXzt4Y9IccTx00sW6/h4FdPSk+K5cFkRZQ0uXttZB6mF4GoAT3d4GtDsD+TSS6wewbRCDeSUOtR+4BgRSRLrP6xTga0RbtOEVNVsZaIszEgcdL3cQCDXpok/R+uv75SxsbKF7563gHSnDk+dCDSQUzHhzld2sbGyhZ9cvIicFEekm6NUxC0tTmdJcTqrt9RQZ1IBYw13DIfmMms4Z7z/witjqtVL5/OFZ/9KRZCIXOwv1XQs8IyIvOBfXigizwIYY94FHgXWYZUesAH3RKjJE1qdf6jklLTBA7nA0Mo6HVo5KhVNLn754nZOmZvLeUsKIt0c5aeBnIp6GytauOM/u7h4eRFnLdKTj1JgDbG8cFkhKYlxPL3DjdcY6GwKz85bynvnxoEVyHnc0F4Tnv0rFUHGmMeNMcXGGIcxZoox5kz/8ipjzDl91vueMWaef0jmp4wxGoEEEehhC/S4DSTNGUdCnE3nyI2CMYbvPrkZY+CHFy3SUU0TiAZyKqq5e7x89eH15KQ4uO38yZVlUKnxlpQQx6VHFFPpdlDR2BmeQK6zGdwtkFHSuyx1ivXsqh///Sulokqghy1viEBORMhLdWiP3Cg8s/EA/9lWy9fOmDPkXEQVXhrIqaj2yxe3s7O2ndsvW0J6ko7nVqq/OVNSWTK7lOrWLpoaasd/hy0V1nPfHrnEDOu5s3n896+Uiip1bV0kxttIcQxdUSsv1aFz5Eao2dXNbU9tZklxutaMm4A0kFNRa82eBv705l4+eXQpJ82JzZTMSg3H6QsL6YlPYce+MCQc6agDYF9XKvevKePOV3bxWpkbxG7VklNKqRGoa+siN9UxrOF+eamJmrVyhH7y7FaaXD387JIlWjNuAtJfREWlNncPX3v4I6ZmJfHtc+ZHujlKTWhJCXHk5uZTX19j1ZYbT64GXL447nuvhvJGF8YYnt9cQ3lngvbIKaVGrLat62Aik6HkpTl0jtwIvL2rnofXVnDDiTNYUJgW6eaoIDSQU1Hp+09v4UBLJ7/6+DKShzHcQqlYV1pcRLK3jdd31o3rfrrb6llXJyTE2/nCybP4r5NnsbAwjY0NQmM4hnYqpaJKXVsXucPMRp2X6qClswd3j3ecWzX5uXu83Pr4RqZlJ/HlU2dHujlqABrIqajz/KYDPPpBBTefMosjSjMj3RylJgVnWg6lzm7W7WuipbNn3Pazt7yCRpPC1UdPJT0pHptNuHxFMV5HOhUHDozbfpVS0amuvWvIjJUBgfU04cnQfr16B2UNLn5yyWIS4+2Rbo4agAZyKqrUtrm59bGNLC5K54t6B0mp4XNmUpgWT5zXxVu7xid7ZIurm8a6aqZMKaQ0uzfzmSPOztSiQjpaGqht6RyXfSulok+Xx0uzq2fYgVxgCKYOrxzcpsoW/vTGXq48qoRVM3Mi3Rw1CA3kVNQwxvDNRzfg6vby648vI14n5So1fIkZJMbbODIX3tvbiKvbE/JdvLFlP3ZfF/NnTjvss9lTi4jDyzvby4e9PVe3h/au0LdTKTU51Ld3A0OXHgjo7ZHTzJUD8Xh9fPNfG8hKTuDWszXHwESnV7oqavz9vf28sr2OW8+ex6y8lEg3R6nJxWkNQz62yEaXx8fbuxpCuvkWVw9b95SRm+ogLTPv8N2nZpOdnMDu/ZV0eYaev7KpsoUfP7OVnzy7lXd2h7atSqnJITBEctg9cmnWetojN7C7X9/D5qpWfnDBQi3bNAloIKeiwr76Dn70762cMDuHTx87LdLNUWry8Qdy2bYOFhSk8vbuhpAmBHhrdz1OTwuF6U5Iygqy/wxyUhKwdbWwpap10G01dnTzj/fKWB5fxim29Ty7voyKJlfI2qqUmhxGGshlJzuwic6RG8j26jZ+s3oH5y4p4OzFBZFujhoGDeTUpNfj9fGVh9eTEGfj55ctxWYbupaMUqqfhGSwxUNnEyfPzaOzx8v7+xpDsmm3f1sLMjw44m3gDBLIJWaQmhhPbnwn68ubB93eGzvrSOup4/yeFziF95nWvYv/bNOMl0rFmkBx7+GWH7DbhJwUh9aSC6Lb4+Prj3xEWmI8P7hgYaSbo4ZJAzk16f1m9Q4+3N/Mjy5aRH768E7mSql+RKxeuc5mSrKSmJmbzJs76/F4fWPe9Lr9Tbh7fCzJ8oHdYQWN/TnSEJuN+Rledta20+oOnjmzvcvDB2VNrMxoxxFnI84mHJvZwtYDbdS0Hj7vZVt1Kz99bit3vrKLpo7uMR+LUmriCPSsZackDPs7Vi05nSPX30+e3crGyhZ+fPEisodZzkFFXkgCORH5s4jUisimAT4XEfmdiOwSkQ0ickQo9qvUmzvr+f2ru7nyqBLOX1oY6eYoNbk5M6CzCYCT5uTS6vbwUUXzmDbp8xne3tVAaVYSOfZ2a1ilBOk1t9nAkcasNC/GwIbylqDbe3tXPR6f4Yi0FohPhpy5zIpvIM4mh82Va3P38I/3yomzCXVtXfxzbTnGmDEdj1Jq4qhr6yIrOWFEyc3yUhN1jlw/a/Y0cN/b+7hm1TTOWqRDKieTUPXI3QecNcjnZwOz/Y8bgD+EaL8qhtW1dfHf/1zPrNwUvne+DgNQasycmQcDuVl5KRSkJ/L6jvoxBT/bqtto6Ojm+Fk54GoKPqzy4P4zSDPtFGc6+XB/02Efd3m8rNnTyIKCNNLcVZA5FbJmkOCqZnmh9Z3O7t55fc9tqsbrM3z2uOmcu6SAsgYXGyt7A8SaVjd3vrKL37+6i9ogvXlKqYmttq1r2BkrA3JTHBrI9eHu8XLrYxspyXJyy1lzI90cNUIhCeSMMa8Dg02muBD4m7GsATJEREN+NWo+n+GrD6+nzd3DHZ84AmeCFqtUasycmdDVBl4PIsIJs3Oobetie03bqDf51q56MpLiWViYBp2NwROdBDhSoauN5aUZVLW4qW45NLh6f28TnT1eTpqRBm3VkFEK6cWA4bgpPXR7DR+UWQFgWUMHH+5v5vjZOeSkODiyNJOC9ESe31RNj9dHe5eHv769j8aObhrbu7nv7X10+EsZNLu6eXt3vSZEUGqCq2sbfjHwgLw0Bw3tXXh92jvv9Rm+9a8N7K3v4KcXLyEpIS7STVIjFK45ckVA3+JAFf5lSo3K3a/v4Y2d9dx2wULm5qdGujlKRQdnJmDAbfVaLSnOIN0Zzxs7RlcgvKq5kz31HRw7Ixub1w09roPZMYNypEFXG0uKM7AJvLu3d6ikx+vjzV31zMhJpsThstqZmg9J2QBMietgRk4yr26vpcXVw+MfVpLujOfkubkA2GzCuUsKaHL18Pi6Sv72zj7auzx89rhpfGbVNNrcHh56bz+bKlv43cu7ePqjA9zxn53s8AexFU0udta04dOLP6UmjLq2LnJHOJ8rL9WBz0BDu96o+d8XtvHE+iq+ceZcjp+thb8nowmV7EREbhCRtSKytq6uLtLNURPUB2VN/OLF7Zy7pIArjyqJdHOUih6BIMvdDFgZ3o6flcOe+g7KG0ee3v/NXfU44mwcNS0LXP5BG4P1yCWkQHcHKfE2VkzL5P19jQcTlLy9u4GWzh5Ompt7cPgnzqze7bmaOG9pAZ09Xn72/DZq27q45IgiHHG9vfUzc1M4flYOH5Y3U9XcyRUrSijOTKIkK4mLlheyu66DB9/dT0ZSPJ87YTo5KQ4eWFPGXa/t5s5XdvPnt/bx+1d3UdbQwY6aNj7c3zSsmndKqdAzxlDXPvIeuVx/hstYH165p66de9/YyxUrivnCKbMi3Rw1SuHqQ60E+l5xF/uXHcIYcw9wD8CKFSv0tqc6TIurhy899CGFGYn89JLFSLCkCUqp0UnMsJ47e+enrZiWycvbanhjZz2fOLp02JtqdfewoaKZldOzraHPjX2Cr4E4UgED3e18bO4U1pU18+8NVZw4J5f/bKtlfkEqc6akwl5/UOjMhDiHlfSks5GCdCefO2EGH5Q1sago3VoXoHoT2BMgdw7nLM5naUk6yQlxZCb3Zro7cmoWyY44Wlw9HDE1k3i7jWuPn87jH1ZS2+rm7EX5pCTG8e+PDnDXa3sOfi/FYWdaTjLuHh/Tc5JYWJhOvN1GujMeu5ZCUWrctHZ66Pb4RjW0EgKlC9LHoWUTn7vHy3ce34Qz3s43zpwX6eaoMQhXIPcUcLOI/AM4GmgxxhwI075VlDDG8M1/baCm1c2jN60iLTE+0k1SKro4M6znPoFcYrydo6dn8frOeho7uslKHl6a7zW7G/AZWDXTGvoY6OUjcZALJ4c/8OpqIz09jdMXTOG5TdVsOdBGisPOeUsKe9tni+9dP6m3x29aTjLTcvqUN2itgvf/aL0+/YdIYhrFmUlBdz8vP+2Q98mOOK4+Zuohy2bnpbCnroPUxDgM1hzAmhY3cXYbq7fW8tIWq56dI85GaVYS2SkJFKQ7Kcp0kpfqGFF2PRWdRORy4DZgPrDSGLN2gPW+DHwOEOCPxpjfhKuNk0FduzWHdsSBnH/9sdaSM8bwwLv7OX3+lElV+qjV3cOn7n2Pj8qbuf3SxSP+81MTS0gCORF5CDgZyBGRCuB7QDyAMeYu4FngHGAX4AI+G4r9qthy75t7eX5zNd8+Zx7LSjIi3Rylos/B3q1DM0YeOzOHN3fV88bOOi5cNvT05lZ3D2/vbmBhYRo5gfkr7hYQq8TAgAKBWXc7ACfOyaU408nuug6Om5XdOxHf1Wj1xgV65J2Z0F4TfJu1W3tf1++A4hVDtn8wqYnxLO1z/pmZm3LwdWNHN2UNHXh9hsrmTsoaXOxvdLHGYwWZNoEpaYnkpyWSm+ZgSmoi+emJZCbF6+iC2LIJuAS4e6AVRGQRVhC3EugGnheRfxtjdoWniRNfYGjkyIdWWuuPNZnR2rIm/ueJTexv6OA75y4Y07bC6dcv7WBjRTN3XX2ElhqIAiEJ5IwxVw3xuQG+EIp9qdj0zu4GfvrcNs5elM/nTpgR6eYoFb36lCAISHfGs7wkk7X7mjhmRjZT0ga/+/zi5hq8PsNZC/N7F7pbrEDNNkiPVJ8euYAZuSnM6BMsAVb2y75JU5KyrIDNmMNr1DXsgpQp0NUekkBuMFnJCQd7LAN7McbQ2NFNZXMnB1rcBxPAfFjefPB7jjgbeWkO8lITyUt1kJfmIDfFQWZSAjYdnhl1jDFbgaGC9/nAu8YYl3/d17CCv/8d9wZOEoFALC91ZL1hjjg7GUnxY54j9+jaCgDe2Dm6ZFCRUN3i5sF393P5kSUaxEUJzTOqJrwDLZ3c/Pd1TMtO4ueXL9U710qNJ2dGb2KSPs5clM+26lb++X45/3XyTOIGGCJY2dzJuv1NnDArh+y+2eQ6mwcfVglWshM4JJALytUI+Yv7tDkTfD1WT56jXxbbtmrImm4FksF67VoqweOG7JmD73OURITsFAfZKQ6WFPcud/d4qW3torrVTXWrm9pWNztq2g6WTwCIswk5KQ5yUx3kpCSQnZJAZlIC2ckO0pxxei6MbpuAH4tINtCJNappoCGYN2DV6KW0dPjzWCe7ulH2yIE1vNKaIzc6rm4P/95QhTPezrbqtlGVQQi3V7bX8vPnt+PzGW7+mCY3iRYayKkJrcvj5cYH1tHl8XH3p1aQ4tC/sir0ROQs4LeAHfiTMeZnEW5S5DgzoWH3YYtTHHFcckQxf3unjJe31XJm3942P2MMz2yoIjnBzinz8g790N0CKbmD7zveCba4wQM5rz9g65s0JZCkJdDrF+DzWXPzkrLB7oDqDYdv641fgPHBSd+CtPDdoU6Mt1OanURp9qHz9VzdHuraunof7V1UNXeyuaqFvpUP4mxCZnICWUnxZKU4yPb3BmYmJZDmjMMZb9dAL4JEZDVw+D8S+I4x5smhvm+M2SoitwMvAh3AeiBoitRYTRRX19ZFQpyNtMSRXxfkpo6tKPiaPQ10dHv51tnz+Nlz23hrVz0XLZ+4VbW2V7dx7X3vk5+WyK8+voySrODzhNXko1fFakL7/tNb+KjcGss9Ky9l6C8oNUIiYgfuBE7HqnH5vog8ZYzZEtmWRYgzAzyd0OOG+EOHLM0vSGPF1Exe21HHvPxUpmb3JhUxxvDMxgPsrXdx0bJCEuPth27X3QI5swfft4jVKzdYINfZbD337d0LvHa3+AuEB9ZtsoK0pGyIS7QCwG4XJPgvYmq3WJ8DVH4Aaef1frf1ACTngj28/00mJcQxNTvukD9bsAr3Nru6rQLm/keD/3lfg4suj++Q9ePtQlpiPGnOOP9zPGmJ8aQmxpHssJOUEEdyQhzOBDsJcZqAJdSMMaeFYBv3AvcCiMhPsM5Pyq+urYu8VMeobljkpSby3t7DRx4MV3ljJwCXLC/iD6/u5t29jRM6kLvzlV0kxdt55ksnDDthlZocNJBTE9Y/3tvP39/dz00nz9Sx3Go8rQR2GWP2APiz614IxGgg55971tkE8Yf/uzt3SQF76tt5eG051x8/g4ykeIyBxz+sZG1ZE6tmZrNyer8SA54uKzgcamglWD1qgwVyXa3+dmb0LusbyPXl8hcUT8qGeH/w1lEHCf5MlC0VVgKW1EJo7C0pQM0WeO9uyJ0Px9zYu9zbA/bIZMu123qHaPZnjKGj20tTRzdNrm5aOz20dPbQ5u6h1d1DeZOL1ioPngGKmcfbhaSEOJIS7Djj7TjibSTYbf5nO444GwlxtoPPcTYbdpscfMTZBLuA3Wb1FNpsgk0EwYrNBQHxvzYGOeQzi8jh0xv7HGDQxcLQF/A2mxAXPzkvXEUkzxhTKyKlWPPjjol0myaS2jEMZ8xLdVDX1oXx/30cqarmThLibOSmOpiekzyqOpvh0OP1sa++g39vqOJzJ87QIC4KaSCnJqTXdtTxnSc2ceKcXL5+xtxIN0dFtyKgvM/7CqwyKQPyeDw0NDSMa6Mixm0Htw0O7Iee4P/pnzEzmYfXlvPzpxpxxttxJthp7Ojm2BnZHFuUQGNjvzvdHQ3WNrviYKg/N48TXC0Dr1dTaW3L5etdx+e12l1bDSl9vndgv7Wu2wbeeOu5ah/4Uno/l1xIKIID66GuzkrGsuV1a93y7VC4E5KzYMeLULUOFl0GObOsY2opt+bq2ezQ3WEFi6MZzmiM9fB0gfGCz2PN+fP2gNcDxuN/9lmvfT5rHeO1jt3/SDJeknxeinyBdb2AFxweTLyXHo+Hnu5uerwePD0evF4PPZ4evF6Dx+ulx+fD5/Hi8/nwGWM9+1+7MbiNzx86GQT/awOCb7Cjiyhb0XIWn/qJSDfjECJyMfB/QC7wjIisN8acKSKFWEO7z/Gv+i//HLke4AvGmObItHhiqmvrYmr26IYI5qY66Pb6aOnsISNp5MFNRXMnRRlORITiTCebq1pH1Y7x9EFZE5++9106ur2kO+M1UVyU0kBOTTibKlv4rwc+YM6UVO78xHItqqsmhJhJKJDoLw/QNfCFSXFmEp8+dhrljS5qW7to6OjmlHl5HFGaGfwL3f4etv6JSIJJSB64lAAcLE1wMDEKWIFUQtLhPXnu5t6SB8YLyKGJXDrqICUf0oqg4n3oqLUyXLaUQ3optOyHxp0Qvxgq11rB0d5XIXMarH/A2p+rAVLzYfMTkF4CpUfDjhcgOceax9e4C9KKrfV8Hqvdrnor6HO3WV1SPSG6my82a46h2K2ANPBa7IjNToLNTkKCHWwJ1p+Z2Pyf23oD0IOvxXqNFbb1GIPXCF4f+IwVuhkjeAGvAYMNr8+HMYH3gV40weDvVPP3rBn8n/tXMcba90jCQTOM3jiAlJzCEWw1PIwxjwOPB1lehZXUJPD+hHC2a7Kpa+9ixbQBzjlDyPNn3q1t6xpVIFfV3ElhhrWNokwnL26uweczEyrL7G9W7yA+zsbnji7lzIX5vaVgVFTRQE5NKBVNLj573/ukO+O577NHkapFv9X4qwRK+rwv9i87RP+EAtnZ2eFpXbj5MsEJJHTDIMeYnQ1zpw748aE690KiD/JLIHWIP7fsbGjeAFlZwXu3anogyQb5xYd+npkBcZ2HtnlfB2SmQ64/yUpmOtjbrHU8XWAaoWAlFC+EPY8BzZCYB/Z2mHsu7O2y2t6TBQ4PlBwD5Wug6j8gLZCRDrVvQS2Q7ICuMti5zwqGu8rBtcsKEls3WyUSnFlW8Fc8C3o6wem/Q56QYgVN8YlW8GWLt4Zw2uPBnmAts8f7PwsEavbe97Y4f2A2cS4ilRpP3R4fjR3dIy49ENC3KPicKcO4wdRPVXMnJ82xzivFmUl0e33UtXcNWZolXMobXbyxs54vnzqbr5w+J9LNUeNIAzk1YbS4erjmL+/j7vHy4E2rJswJUUW994HZIjIdK4C7EphYY7HCyebvwepXS25M3M3W83DmyCWkWL1nPS6rd+6wbbVa7esftCSmH142obPRmh8XkJxr9cKBv9fPWJkqk7Kt/TaXWUERQNYM6KiHPa9YQV9SNiy6xBqCWb7G6n07+vPwyk+tYZAnfBWay6F+O8y/wDoGdwtkTLV64mxxGmgpFSINHaMvPQC9gVxd+8hLEHR7fNS2dVGY4QSgONN6rmhyRfS6xRjDz57fRmtnDzYRbAJXHFUy9BfVpKaBnJoQujxebrh/LfsbXPz12pWjukOm1GgYYzwicjPwAlb5gT8bYzZHuFmRFaQo+Jh0NltZI+OHcZHjCAztbA8eyHW1Bg8IEzOgce+hy1wNMGVR7/vkXGuemzFWfTmA1AIrwMoohaYyQCA+2Rpimb8Idr9sBXizToc4B8w+3ZovN/8Ca6joKbda20tMs4ZYlhzVu79A4pgIJUhRKlqNpYYc9Bla2TryEgTVLW6MoTeQywgEcp0cOdxRCuPg1R113P1ab9KmsxbmU+Rvm4peGsipiPP5DF9/ZAPv7m3kt1cu49iZUTpkTU1YxphngWcj3Y4Jw5lpBS+h4m4ZXm8c9M6j62qD1CnBt5UapDxXYjr0dPRmlvR0Wdvo2yOXkmv19HV3QNsBq5csKcf6LGMq1G615uBlTfcHd9Os77tboWSltd6s02Dmqb29a8OZ96eUCqlAIJc3ykAu2Z+hdTS15CqarTmtgQCuKLM3kBtKU0c3if4EUaH2zIYDpCbGcecnjuDFLdXcfMoQ5V5UVNBATkXc7S9s4+mPqvjmWfO4cNnErcOiVMxwZlrFs40JzXDAUQVyAyRbcbdCTpA5H4HtdzZbAVtgmGVSn1IIyf65ch21Vp24lHxrKClA5lTAWIFc7jxrmc0GJ3zNms+WnNO7HR0iqVRE1Y6xR05EyEsbXVHwqmZrOGagRy4pIY7s5IQhA7ker48VP17N8pIMHr1p1cgbfVg7Ovn6Ix9x4bJCLl5ezIubqzl9wRROnJPLif75eyr6aSCnIupv7+zj7tf2cPUxpdx4kqbGVWpCcGZa87q62nqzWI7FcIqBBzj82SgD2Sn78nQPXI8uUFfO3eIP5PxlCJxBArn2Omirguw+bcqcbiUcMQYKlvYuT0gOPsRTKRUxgR65sWRizEt1UNs68jlylf6ArSCjd6h4UaaTyubBA7mXttTg9RnWljWxs6aN2WOcQnLXa7t5e3cDb+9uoKWzh1a3h/OWaM3dWGOLdANU7HphczW3PbWZ0+bncdv5C0dVlFMpNQ4CQVEo5sn5fP55bRnDWz8+GZDgRcEDvXRB58j1KwreGeiR6zO0MikH7A6o2eRPRNInEUB8Ipz6XTjl26EJXpVS46aurYvMpHgS4kZ/GZuXmngwIByJ6tZOclIScMT1Do8sSE+kumXwQO65TdUEqhO8sr026Do+n6Gz23vY8to2N//zxCZ21faeF9/YWc/SkgxSHHH85Nlt5KQkcMJs7YmLNRrIqYh4eWsNN/99HUuKM/jdVcuJs+tfRaUmjECSjlAEct1tVv214Q6ttNmsXrmuID1ygSBtoGQn0Jsh09Vgpe7vO4fNZrOCt+oN1vuMfpkJnJmQkje8diqlIqa2zT3qYZUBuamjG1pZ19Z9WE9gQbqTAy2D9+7tqm3npDm5TElzsK06yI0q4LanN7PyJ6vZ33BobclfvLCd+9eU8ZV/fgRYGTL31ndw4dJCvnbGHETghhNnEK/XUjFHh1aqsHtlWy03PbCO+QVp/PXalSQl6F9DpSaUUAZygwVfA3GkBZ8jF9iWI0iPWXyi1dsWWKej3uqN69/Tnz0LGnZZr9OLh98mpdSEUdfWNeZALi/NQXuXB1e3Z0TXIXXth+97SloibW4PHV0ekh2Hb8vnM+ytb2fVzGx8BrYdODyQq2vr4m/vWEmm/vzWXm67YOHB767eavXgbaxsYVdtO+vKrHPz8bNzmDMllY8fVaLXUjFKQ3cVVqu31PD5+z9gTn4K9197NOlOTcut1IQTn2T1ZgV6t8ZiNIFcQsoQQysHGProzOgNPl39asgFzDjZqhG39BNaFkCpSaquvYvcMcyPAw4WEx9pCYL6tsP3XZBubat6gDl3B1rduHt8zMhNZl5+Krtq2+nx+g5ZZ315MwBJCXbe3l1/cPmO2jYaO7r5+hlWkqeXt9bwxq56clMdzM5L8X9Hg7hYpYGcCpvnN1Vz04MfMK8glQevO4b0JL2IUmpCEgldLbnOZut5RD1yqcEDOXerVTIgPin49wIFv30+KzNlsGGS8U447stQevTw26OUmjCMMdS1dR2sBTdavUXBhx/IGWOob+8ip1+PXH4gkBtgeOWeOmuo+MzcFGZPSaXb62N/46HDJzdWNGO3CdcfP50dNe00u7oB2FRp3cA6a1E+CwrSeGFzNW/vquf4WTmaW0BpIKfC49mNB7j57+tYWJjO/dcdrUGcUhNdYkbohlaKLfhwyIE4UoJnrXS3WNsZ6OIlZQq01/qDOY9VXkApFVXaujy4e3xj75FLs74/kh65ti4PXZ7D9x3okRtonty++g4AZuQkMz3HuhHVfx7c5qpWZuWmsHK6NZIgEMBtrmrBGW9nek4K5yzOZ93+Zho6uvnYPJ3PqzSQU2Hwz/f388WHPmRZSQb3X7dSh1MqNRmEqkfO3WL1sNlG8N+NIxW83VZR7766Wgfv2UvNB+O1slJC8ILiSqlJrW6MNeQCAsFYbdvwSxDUB8oepCYcsnxKWqBHLnjmyspmNwl2GzkpDqZmW+VM9jV0HLLO3voOZuYls6jIuum1sdIalr6psoUFhWnYbcJVK0spynCyoCCNMxbq+U2FKJATkbNEZLuI7BKRbwX5/BoRqROR9f7H9aHYr5rYfD7Dz57bxjf/tZFVM7O579qVpCZqEKfUpODMtIY3ej1j285IioEHJASKgvcbXuluHbw0QIr/wubA+kPfK6WiRqCOW2A442hlJiUQb5cB57UFczCITDl034nxdjKT4gfskatq7iQ/PRGbTchOTiDFEUdZnx45j3+o5bTsZDKSEijJcrKpsgWvz7C5qpVFhdZ5LzvFwX++fhJP3XzcIeUPVOwa8+xIEbEDdwKnAxXA+yLylDFmS79V/2mMuXms+1OTQ2e3l68+vJ7nNlXzyaNL+f4FC7XEgFKTSSBzpbsFkoMkDRkud8vIU/oHSgZ0tUJyzqHbGqyweGq+NYyzeT+kFlrz4ZRSUSXQkzU9J3lM27HZxCob0DyCHrl2a95a/x45gPx054Bz5A60dB4cfikiTM1OOqRHrqKpE4/PMM1/TIuL0tlQ2cze+g5c3V4WFfXeDNMATvUViivrlcAuY8weY0w38A/gwhBsV01StW1urrznHZ7fXM3/O3c+P7pokQZxSk02oSpBMJoeuYOBXJ95cp5u8HQOvq04BxQstV4XrxjZPpVSk8Le+g6SEuwHk5WMRWFGIlXNgxfy7qvOPwwz2Py8gvTEAXv3qprdFGX03lialpN8SI/c3n7B6ZLiDMobO3nVXzi8byCnVF+huLouAsr7vK/wL+vvUhHZICKPikhJCParJqBt1a1cfOfb7Khp5+6rj+T6E2ZoViWlJiNnhvU8lkDO0+UPvjJG9j1HkKGVB0sPDHFBs+hSmH4STD1uZPtUSk0K++o7mJqdHJJri8IM54gCufr2buw2ITMpWI9cYtAeOa/PUN3qpiCjdzjmtOwkyhtdePwlCPbWHRrIHTvDGgXxm9U7SUqwM8tfZkCp/sLVTfI0MM0YswR4CfhrsJVE5AYRWSsia+vq6sLUNBUqr26v5bI/vEOP18fDnz+WMxZqxjilJq1Q9MiNpoYcWHXk4NBA7mAZg4zBv+tIhUWXWAXClZqAROTnIrLNf3P7cRHJGGC9QfMPxKp9Da6DmR/HqjDdSXWr+2BANZS6ti6ykxOw2Q4PIgvSEmno6Mbd4z3sO16foSC9t0duanYyHp+h0h9E7mvoINURR3ayFSAuKkonJ8UqWH7q/CnE66gmNYBQ/M2oBPr2sBX7lx1kjGkwxgTSj/0JODLYhowx9xhjVhhjVuTm5oagaSocjDHcv6aM6/66lpKsJJ68+TgWF+swAKUmNXu8FVBFIpCzx1kBWd99u/wFcoMV+VZqcnkJWOS/ub0DuLX/Cn3yD5wNLACuEpEFYW3lBOTx+ihvdI15flxAYYYTn4HatuGVIKhv7yJngLIHU9KDFxiv8meyLDykRy6QudIaXrm3voNpOb29jHab8OOLF7FyWhZfPnWQecEq5oUikHsfmC0i00UkAbgSeKrvCiJS0OftBcDWEOxXTQAVTS6u/+ta/ueJTZw0J5dHbjz2kLtOSqlJzJkB7ubRfz8QyAWGaY5EUk5v8AbgarASmYxmW0pNIMaYF40xgXSwa7BugPen+QeCOJgUJDtUgZwVXA13eGWwYuABvbXkDt1WYNuFfefIZVs9ioH6cnvqOg4LTs9cmM/DNx6rwyrVoMactdIY4xGRm4EXADvwZ2PMZhH5AbDWGPMU8CURuQDwAI3ANWPdr4osj9fHX97ax69e2gHAt8+Zx7XHTdekJkpFE2cmdNQPvd5AAkHgSHvkwMpWWb+z972rwRpWadOMbSqqXAv8M8jyYPkHjg5Liyaw/klBxiqQgKSyuZPhpEeqa+tiVl5q0M8CgVz/hCeBrJh9b3LnpjpITYxjR00bbe4eKps7+cTRpaM4AhXrxhzIARhjngWe7bfsu31e30qQoQNqcvpwfxPffnwTWw+0cuq8PL5/4UKKM0MzXl0pNYEkZhwaTI1UZzPEJVrZJEcqKQfca606dvY4K5DrW4pAqQlMRFYDwSaKf8cY86R/ne9g3eB+cIz7ugG4AaC0NLqDgUAP1rQQBXIF/kCuahglCHw+Q11714CFyPP9gVr/hCdVLZ0kJ9hJS+y95BYRFhSksbmqlR011lzgOVOCB4hKDSYkgZyKDR6vj9+s3smdr+5iSmoid119BGcuzNeslEpFK2cmeNzQ0zm6mmzulpFnrAxIzgGMFcClTrGe8xaObltKhZkx5rTBPheRa4DzgFONMSbIKkPmH+izr3uAewBWrFgRbFtRY1/9oUlBxirFEUe6M35YQyubXN30eA35acEDuRRHHCmOuMOKglc1d1KQ4TzsWmlhYTp/f6+MzVVWRt55+RrIqZHTQE4NS0WTiy//Yz0flDVx+ZHFfPf8BaQmxke6WUqp8dQ3c+WoA7lRJj5KyrKeXQ1WO7raNNGJigoichZwC3CSMcY1wGoH8w9gBXBXAp8IUxMnrL0NrkOSgoTCcEsQBIZMTkkbOCNusBIEB1rch8yPC1hcnIb7LR9/eHU3U9IcFGdqfgE1cjqhSQ3puY0HOOe3b7C9uo3fXrmMn1++VIM4pWLBwUCueXTfH1Mg5x9G6aq3gjmAZA3kVFS4A0gFXhKR9SJyF4CIFIrIs2DlHwAC+Qe2Ag8bYzZHqsETxd769pANqwwoTE88WAZgMIFslIHslMEUpCcGSXbipjDId06dPwW7TTjQ4ubU+VN0dJMaFe2RUwNqdffwg6e38OgHFSwtTud3Vy1naogyRSmlJoGx1JLz+awi3qMN5BypYHdYyVYCwzOdWaPbllITiDFm1gDLq4Bz+rw/LP9ALOv2+Khs6uTiZUUh3W5hhpO1ZUOf42qG0SNXnJnEc5UHDr7v8nipb+8Kms07LTGez50wg7d21XPtcdNH0XKlNJBTQXR0eXjw3TJ+/+puWjt7+NLHZvHFU2drQUqlYo0jzUr5P5pArrsNjG/0gZyINZTS1WBtx54A6cGytCulYsH+Rhc+E7pEJwGFGU5aOnto7/KQ4hj4sjgwtDJvgGQnADNzk2l29dDY0U1WcgL7/XXipmYHTwj3rbPnjaHlSmkgp/zcPV5e2VbLvzcc4OVtNbh7fJw4J5dbzpzLoiIt7q1UTLLZrN6w0QRyY6khF5BWCLVbrCAuZ45VpFwpFZNCnbEyIFBL7kBzJ7MHyRxZ09pFTkrCoDe1Z+ZaNd/21LWTlZzF7rrQlktQqj8N5GKYMYYPypr4+7v7eWFzNR3dXrKTE7jsyGIuXl7MkVMzI91EpVSkOTPHFsiNtkcOYMZJULkWelwwa9AkgEqpKLcvUEMuxFM8+taSGzyQc5OXOvCwSoAZuVbb9tR1sGJaFnv9wWdguVKhpoFcDPL6DE+ur+SPb+xl64FWUhPjOH9pIectKeSYGVla1Fsp1cuZAU37Rv69QIKUsQRyGaUwZRHUbIK8BaPfjlJq0ttb30G6M57MEJUeCCgcZi25mlY3+YMkOgFrjlyC3cbuunbA6pmzin/raAI1PjSQiyHGGF7ZXsvtz21ne00b8/JT+cnFi7loeSFJCfpXQSkVhDMTqtZbyUtsI7jJ426x5tcljLE20uLLoegIzVipVIzb19AR8mGVYM15swlDliCoae1iSfHgN6bsNmFqdtLBIZV76jt0WKUaV3r1HiM+KGvk9ue3897eRqZlJ3HnJ47gnMVazFspNYTEDDBeKwPlSOa7uVusZCkjCf6CcWZA0ZFj24ZSatLbV+9i5fTQZ66Ns9soSHeyv3Ggkn7Q4/XR0NE15NBKgDn5qawra8Lj9bGzpo1zlxSGsrlKHUIDuUjwdEPNRih/H6rWWRnZHKnWRY8jFaYeB6XHWFnbxmhbdSu/eGEHq7fWkJvq4IcXLuTKlaWagVIpNTyBEgTu5pEFcp1NY0t0opRSfu4eL1UtnUwbpxJIs/JS2FnbPuDndW1dGMOQQysBTpydwzMbDvDw2gpa3R5WzdTRBGr8aCAXLg27YcPDsPc1qPoQPP6x2KmFEOeArjbr4bUKTpI1E5ZfDUuvgrSCEe2qs9vLK9tr+fu7+3lzVz2pjji+ceZcPnvcNB1CqZQamb615DKnDf97rgbI0tpISqmx29/owhiYlhM8jf9YzZmSwjt7GvD6DHbb4TfRe2vIDVx6IODkuXkAfPvxjdgETpydG9rGKtWHXtWPp44G2PwYfPQPK/Oa2KDwCFhxHZQcBcUrIb1fYcuuNtj6NHz4ALz8fXj1Z3DyN2HVl4Km3u7yeKls6qSs0cXmyhbW7W/mnd0NdPZ4yU9L5BtnzuWTR5eSkRTaycFKRQMRuRy4DZgPrDTGrI1siyagg4Fc8/C/4/NZPXhJeidaKTV2geyP4zXfbPaUVLo9PsoaOpjhLyHQV83BGnJD98hNSUtkUVEamypbOXJqJulJmuhEjR8N5EKtxw07nocN/4SdL4LPY2VdO/2H1qT9IXrX2nFSmX8+VatOo33aduZt+jWzX/4BlW89xJ+yvsZm3zRa3T20dFoPV7f3kO/PykvhsiOLOXtRPiunawZKpYawCbgEuDvSDZmw4p1gd4ysBEFnkzVkXAM5pVQIjFcNuYA5/rIDO2vbgwZye+sHL+zd328+voz39zVxzAw9B6rxpYFcKPh8UL7G6nnb/AR0tUBKPhxzEyy5EvIXBfmKYU99Ox+Vt7CtupUdNe3srGmjqqV/+tvruDBhKd91/5nvVH2BB1Ov482cK1hUlE66M550ZzyFGU6mZicxJy9V7/woNQLGmK2AJv0ZjMjIa8m5GqxnZ+gTEyilYs/e+g6ykxNIG6c0/rPzrOBtZ00bZy7MP+zz3XXt5I2gjMCsvFRm5Y0xY69Sw6CB3FjU77SCt40PQ/N+iE+G+efD0o/D9JPAZj+4am2rmw/Lm/movJmPKprZUN5CW5cHgIQ4GzNzUzhqehZzpqRSkpVEYXoiBRlO8lIdxNvPhc4vwZM385lt9/CZkho4/w5ITBt92709VpHdhJRD2qmUUocZaSDX2Wg9a4+cUioE9taPT+mBgGRHHEUZTnbUBE94sruunZlBeuqUijQN5Eai2wXl70LZW7BrtZW0RGww4xT42P/AvHMhIZk2dw+b9loB2/r91vMBf09bnE2YV5DKBcsKWVaSwdKSDGbkJA89BNKZCR9/AN7+Hay+DWq2wMfvh7z5IzuGht3w+i9g2zNWz2F8shV4nvRNSD38LpRSk52IrAaC/eX+jjHmyRFs5wbgBoDS0tIQtW6ScGZCS/nw13c1WOfGwPw6pVRYtHb28PymA32WHDraoO/gg7TEeI6ZkTUpRiTsa+jg+FnjmzRkzpQUdtS0HbbcGMOeug7OWzKyxHNKhYMGcoPp7oDy92Dfm9aj8gPw9YDYoXA55owfUVVyHptanWw90MrWj7ay5UAr5Y29RSWnZidx1LQslpZksKwkg4WFaSTGj7IHTASO+7KVMOXRa+GPH4PzfwtLrhj6u54ueO1/4a3fWlkyF1wEefOsgPDDB2HLU3Dpn2DmKaNrm1ITlDHmtBBt5x7gHoAVK1aYUGxz0nBmQHe71ZMfJOnSYVwNVv25sdaQU0qNSFmjixsfWDfs9W88aSbfOnveOLZo7FzdHmpau5g+ThkrA2ZPSeWtXQ14vL5Dbq43dHTT0tmjPXJqQtJAri93ixWsHRK4eTBipyN7MZUzPs0WxxLe98xlZ4th24tttLm3AFaMNT07mSVFGVx5VCkLCtNYWpxBVvI4ZIucfgLc+AY8cg089jnY/hyc8aPDM2AGlL0Dz3wVardYc/ZO/wGkTun9/Pj/trb1wKVw3q/gyGuCb8fTBXXbrVp3mlZcqdjRN3NlyjDuirsadFilUhEwOy+FR758AgCm3+0mw6ELHlizn7te201OSgLXnzAjXE0csX3+RCPjObQSrD+7bq+P/Y2uQxKe7KmzEq3MyB3f/Ss1GiEJ5ETkLOC3gB34kzHmZ/0+dwB/A44EGoCPG2P2hWLfw2IMdLVCaxW0VEJrJZ7mCnqaKjAtlUhbFfEdB4jrscZGe7GzJ2EO79kv5OXuObzrmUVHhRMqrM1NSeugNCuJC5cVMr8gjfkFaczLTw1vjbbUfPjM0/DGL+HNX1uZMldca83Ry18M3m7Y/y6sf8AqZ5BWBJ94GOacefi2cufCtS/Ao5+Fp79szf079btWzx1YAe5798A7v++d+1K8Ei68E3LnhO+YlQoxEbkY+D8gF3hGRNYbY4L8I4lxiRnWc2fT8AO5vIXj2iSl1OES4+3MLxje/PkfXbSIls5ufvTMVrKSE7jkiOJxbt3o7GvwZ6wcp2LgAYHMlTtqDs1cubvOujbUHjk1EY058hARO3AncDpWqPO+iDxljNnSZ7XrgCZjzCwRuRK4Hfj4YNt1d3Wzdfe+YbXB5/PQ01KDp6kC01qBve0A8R3VON3VpHTVkN5Th9N0HvIdmxFaSeeAyaLaZHPAzKDaZLHNlLLLsZDstGxKMpOYneXklMwkSjKdlGQlUZThHP3QyFCzx8PJ34KlV8JL34N374Z37jh0HWcmnHiL1euWMMhJMDENrvonPP9Naxvbn7WCQnerVQvP3QKzz7Tm07XVwBu/gHtOgjN/DEd+tnfgfbfLSvySXmT13Ck1gRljHgcej3Q7JrxAj5y7eeh1vT1WPUztkVNqQrPbhF9/fBnNrve55dENZCYlcMq8vEg36zB7x7n0QMAsf+bKHTVtnLWod1r17tp2HHE2ijKc47p/pUYjFF1IK4Fdxpg9ACLyD+BCoG8gdyFW0V2AR4E7RESM6d/x3yuxYTPz7186qgZ5jVBLJnWSTYW9iBbHctodebid+XQn5eNLLcSWVkBKchLpznjSEuNY5ownzRk/ovSyE0bmNLjir1awtfMlK5CKc0D2bJhxMsQNc3inPQ7O/aUVsL3xS3jrd1bwN/t0OO6/oXBZ77oLL4YnboR/fwU2/gumroID62Hv6+BxQ5wTVl4PJ34DEtNDfshKqTByZljPw8lcGSg9kKSlB5Sa6Bxxdu7+1JFc9cc13PTgBzx4/TEcOXViJSnaV99BbqqDFMf4jnpKdsQxLz+VN3fV86VTZx9cvqe+g+k5ydhsEz8pjIo9ofhXUQT0TWdWARw90DrGGI+ItADZQP1AG3UnFbBt+ZeH2QTBlpJHfFYJidnFpGQXMcXppCDW/tElpsPiy8a+nTlnWI9AnB0so1VaAVz9OKy91xraWfamFTge8WkoPgp2vQxv3wEbHoZjb7buzjfvt4LF2WdAweiCdKVUBNjjrVIlnc1Dr+sKlB7QQE6pYETk58D5QDewG/isMaY5yHp/Bs4Dao0xhxekDZHUxHj+cs1KLr/rba69730evfFYZk+ZOCNq9jVYgVQ4nLEwnzv+s5P69i5yUhz4fIaNlS1a2FtNWBMq2Un/9N7zLvxGhFsU44ZKSWyzwcrPWY/+2eyWXAHH3AjPfB1e+p/ABgED//kRzDnb+jxjKjjSrKFbmuFOqYlruLXkDvbI6YWPUgN4CbjVf2P7duBW4JtB1rsPuAMrx8C4yk11cP91R3PJH97m039+j0dvWjVhhhLurXdxapiGfJ65cAq/e3knq7fUcOXKUt7f10hdWxenL5gy9JeVioBQBHKVQEmf98X+ZcHWqRCROCAdK+nJIWI6vfdkFywleeFy+NzL0FZtFR9PK7Ke37/XGra547neddOK4KjrYdGlkF6iQZ1SE40zE9prh17P1QC2eOsGjVLqMMaYF/u8XQMEHUpjjHldRKaFpVFASVYSf/3sSj5+9zt8+t53eeTGVeOTeXsEdtW2Ud/exYLC8JxPFhSkUZLl5IXN1Vy5spSnN1SRGG8LWyCp1EiF4mr5fWC2iEwXkQTgSuCpfus8BXzG//oy4D+DzY9TUSY1H7JmWPP2nJlw4tfhq5vhk/+Ci/4AZ90O2TPh5e/Db5fAj3LhruPhzd9YmUaVUpGXlGUFaT7f4Ou5Gqx1J0GRYaUmgGuB54ZcK0wWFKbxp8+soLypk2vvex9Xtyei7XlsXSV2m3DO4vAU4xYRzlyQz1u7Gmh2dfPcxmpOnT+F5HGen6fUaI35b6Z/aMDNwAtY5Qf+bIzZLCI/ANYaY54C7gXuF5FdQCNWsKdimSMVZvep03zMjVaNun1vQku59bz6e9YjcxoUrbCyaM4+AxLGtyioUiqItCLw9UBHrXVzZiCuRh1WqWKeiKwGgv1D+Y4x5kn/Ot8BPMCDY9zXIdNSxuroGdnccdVybnzgA258YB1/+vQKEuLCP0rG5zM8ub6KE2bnkJvqCNt+z1qUz5/e3Mtpv3qdho5uzl9SGLZ9KzVSIbnFYIx5Fni237Lv9nntBi4Pxb5UFMudaz0CGnZbNfCq1sHe12DTo2B3QPEKmHUazDtP69gpFS7p/hpTLZWDB3KdjZA1PTxtUmqCMsacNtjnInINViKTU8c6Qmk8pqWcsTCfn1y8mG89tpFvPPoRv75iWdizNr67t5HK5k5uOWvu0CuH0JFTM/nl5Ut5cUs1HV1eTp47jNqZSkWI9hWriSt7plX/DsDntXrpdr4I+96whmG+/H0rU+bcs61HydFgmyA1/pSKNilTwBZn9ZgXHxl8nW6XNQ/WObHSlys1kYjIWcAtwEnGGFek2zOQK1eW0tDRzc9f2E5WcgLfPW8BEsYh049/WEGKI44zFgxy42gciAiXHlnMpUdOzALpSvWlgZyaHGx2mHGS9QBoqYDtz8G2Z2DNH+Dt31nDueaeDfPOt+rnxSdGtMlKRRWbHVILoLV/Lqs+XP6KMsk54WmTUpPTHYADeMkfGK0xxtwoIoXAn4wx5wCIyEPAyUCOiFQA3zPG3BvOhv7XyTOpb+/iL2/tIyfFwRdOmRWW/bp7vDy3sZqzFuXjTNAbtEoNRAM5NTmlF/eWPnC3wq7VVlC35Sn48AGr5tXs063hl7PPgETNoKfUmKWXwIH1Vo3JYHfmm/0lRdP0TrZSAzHGBI2GjDFVwDl93l8VtkYNQET4n3MX0OjvmctJSeDjR419Ht5QVm+toa3LwyXLi8Z9X0pNZhrIqckvMQ0WXWI9PN2w93XY9jRsexY2Pw72BKuHbv75MPcc7S1QarTSi2H/21Y9uWAFv5vLrJsoWgxcqahhswk/v2wpTa4ebn1sIxlJCZy5cHyHOz6+rpKC9EQtxK3UELRYl4oucQlWNszzfwtf2waffR6O+hzUbYOnvgi/mA33nQdr7rKGZyqlhu9gwpPy4J8374eMUi09oFSUSYiz8YdPHsHi4gy++NCHvLvnsFLAIdPQ3sVrO+q4cFlR2BOsKDXZaCCnopfNDlOPhbN+Al/eAJ9/A074OnTUw/PfhF8vhHtOgTd+CbXbrOFiSqmBpRWC2KzMlf31uKGtGjKmhr9dSqlxl+yI4y/XHEVxppPr/7aWrQdax2U/T39UhcdnuFiHVSo1JA3kVGwQgYIl8LHvwBfWwM1r4dTvWZ+9/AP4/dHwu2Xw/K2w5zXw9kS0uUpNSPZ4K3tlsN7slgrAWD1ySqmolJWcwP3XHU1yQhyf/vN7lDeGPunm4x9WsqAgjbn5qSHftlLRRgM5FZtyZsMJX4UbXoGvbIFzfwU5c+D9e+FvF8D/zoRHr4ONj0Jnc6Rbq9TEkV4MrUECueYy61kDOaWiWlGGk/uvW0m3x8en7n2X+vaukG17d107H1W0cMkR2hun1HBoIKdUehEcdR188hG4ZQ98/AErMcqeV+Ff18HPZ1rz6t75vVWkXKlYllYE7hYrW2xfzWVWCRBHSmTapZQKm9lTUvnzNUdR3ermmr+8R4srNKNYHl9XiU3ggqWFIdmeUtFOAzml+nKkWEHcRXfC13fAdS/Bqi9CRx28cCv83xHwmyXw1Jdg02PQMX4TvpWakNJLrOf+9eSa9+v8OKViyJFTM/nDJ49k24E2PvbLV3novf14faOfa+7zGR7/sJLjZ+eSl6Z1YJUaDi0/oNRAbHYoWWk9TrsNGvfAztVWT93mx2HdX6318pdY5Q2mnQDFKzT1uopuaf475S0VkDffeu1utUoSTD8pcu1SSoXdKfPyeOILx/H9pzdz62MbeWBNGd87fyErp4/8/8H39zVS2dzJN86cOw4tVSo6aSCn1HBlzYCjb7AeXg9UfWgFdXtehTV/gLd/Z62XM9cfAB5tPWfPBpt2fqsokZBkDaHsm/Ckeb/1nKk9ckrFmkVF6Tz8+WN5esMBfvrsVq64+x3OW1LArefMpyjDOaxteH2Gh97bT1KCnTMWThnnFisVPTSQU2o07HFQcpT1OOkb0N0Bleug/F0ofw+2/Rs+vN9aNz4ZpiyAKYsgfxFMWQx58yAxPbLHoNRopRf3C+TKrLIEacWRa5NSKmJEhAuWFnL6/Cnc9dpu7nptN6u31vD5E2dy40kzcSbYg36vxdXDw2vL+duafZQ3dnL1MaUkJeilqVLDpf9alAqFhGSYfoL1APD5oGEXVLwH1RuhehNsfgw++Evvd5JyrF6+rBmQPRMyp0FqAaTmWw+Hpl5WE1R6CRz4CHo6Id5p9cilFkJcQqRbppSKIGeCna+cPofLVxTzs+e28duXd/LI2nK+dc58zl9SgIhV4HtnTRv3vb2Px9ZV0tnjZeW0LG49ez5nLNDeOKVGQgM5pcaDzQa5c6xHgDFWL0bNJqjbbs25a9wD+96EDf84fBvxyZA6xZ8JMM3qwTv4SIP4JLAnQJwD4hL7vHaALc7qIQk8EP/r/s+Bz2Rsx5s5HeJ1cnrMSPOnBt/+HLgaoH6nNZRYKaWA4swk7vjEEXz62EZue2ozX3roQ+5/Zx9XrCjhyfVVvLmrnoQ4GxctK+Qzq6axsFBHqCg1GhrIKRUuIpBRYj3mnn3oZz2d0FwObQegvcZ6bvM/dzaBu9kavuZusR7e7ogcwoA+/4ZVcF3FhowSQGDva9aNhtJjYOapkW6VUmqCWTk9i6e/eDwPry3n5y9s5xuPbqAgPZFvnDmXq1aWkpWsvfhKjYUGckpNBPHOw3vwBtPjBk8neLrB47YCO4+7973PAxgwPv/D+B++3kffz8cqo2Ts21CThyMVjv+K1Suckhvp1iilJjC7TbhqZSnnLC5gR00by0syiLNrAjClQkEDOaUmo/hEHcqoIkszVCqlRiDdGc9R07Q8j1KhpLdElFJKKaWUUmqSGVMgJyJZIvKSiOz0P2cOsJ5XRNb7H0+NZZ9KKaWUUkopFevG2iP3LeBlY8xs4GX/+2A6jTHL/I8LxrhPpZRSSimllIppYw3kLgT+6n/9V+CiMW5PKaWUUkoppdQQxhrITTHGHPC/rgYGquSYKCJrRWSNiFw0xn0qpZRSSimlVEwTY8zgK4isBvKDfPQd4K/GmIw+6zYZYw6bJyciRcaYShGZAfwHONUYszvIejcAN/jfzgW2A+lAS5/V+r7v+zoHqB/0YIan//5Gu+5AnwVbPtxjnMjHO9DnIz3e/u8n8jHrbzz8ZX3fTzXGTOqc9SJSB5QxOX/HUP297f8+ksc8mr+3Ay2Phn+r+hsPvnygY4zWcxME/x1D9Rv23/5Y1tN/qwMvH82/Vf2NR28i/cYDn5uMMaN+YAVaBf7XBcD2YXznPuCyEezjnoHe93u9dizHMtD+RrvuQJ8FWz6CY5ywxzvcYxvqeCfTMetvPPrfOFoek/F3DNXf24l0zKP5ezuaY54ox6u/8fgeczQ8hvM7huo3nIi/40T5ezvUuuP9b1V/4+j6jYM9xjq08ingM/7XnwGe7L+CiGSKiMP/Ogc4Dtgygn08Pcj7/p+Fwki2Odi6A30WbPlwj3EiH+9An4/0ePu/n8jHrL/x8JeNx3FNBJPxdwzV39v+7yN5zKP5ezvQ8mj4t6q/8eDLY+H8FOu/40Q53qHW1X+rgy/X33gIQw6tHPTLItnAw0ApVlf+FcaYRhFZAdxojLleRFYBdwM+rDl5vzHG3DvqnQ7clrXGmBWh3u5EFWvHC7F3zLF2vNEqFn/HWDvmWDteiM1jjjax+BvG2jHH2vFC7B1z3Fi+bIxpAE4NsnwtcL3/9dvA4rHsZ5juCcM+JpJYO16IvWOOteONVrH4O8baMcfa8UJsHnO0icXfMNaOOdaOF2LsmMfUI6eUUkoppZRSKvzGOkdOKaWUUkoppVSYaSCnlFJKKaWUUpOMBnJKKaWUUkopNclEZSAnIski8lcR+aOIfDLS7QkHEZkhIveKyKORbks4iMhF/t/3nyJyRqTbEw4iMl9E7hKRR0Xkpki3R41OrJ2fYu3cBLF3ftJzU3SItXMTxN75KdbOTRD956dJE8iJyJ9FpFZENvVbfpaIbBeRXSLyLf/iS4BHjTGfAy4Ie2NDZCTHbIzZY4y5LjItDY0RHu8T/t/3RuDjkWhvKIzwmLcaY24ErsCqx6gmiFg7P8XauQli7/yk56boEGvnJoi981OsnZtAz099TZpADrgPOKvvAhGxA3cCZwMLgKtEZAFQDJT7V/OGsY2hdh/DP+ZocB8jP97/5/98srqPERyziFwAPAM8G95mqiHcR2ydn+4jts5NEHvnp/vQc1M0uI/YOjdB7J2f7iO2zk2g56eDJk0gZ4x5HWjst3glsMt/R6Ub+AdwIVCBdUKCSXSM/Y3wmCe9kRyvWG4HnjPGrAt3W0NlpL+xMeYpY8zZQEwMe5ksYu38FGvnJoi985Oem6JDrJ2bIPbOT7F2bgI9P/U1af+h+hXRe/cIrJNQEfAYcKmI/AF4OhING0dBj1lEskXkLmC5iNwamaaNi4F+4y8CpwGXiciNkWjYOBroNz5ZRH4nIncThXeVolCsnZ9i7dwEsXd+0nNTdIi1cxPE3vkp1s5NEKPnp7hIN2A8GGM6gM9Guh3hZIxpwBrzHBOMMb8DfhfpdoSTMeZV4NUIN0ONUaydn2Lt3ASxd37Sc1N0iLVzE8Te+SnWzk0Q/eenyd4jVwmU9Hlf7F8WzWLtmGPteCE2jzkaxdrvGGvHC7F3zLF2vNEqFn/HWDvmWDteiM1jnvSB3PvAbBGZLiIJwJXAUxFu03iLtWOOteOF2DzmaBRrv2OsHS/E3jHH2vFGq1j8HWPtmGPteCE2j3nyBHIi8hDwDjBXRCpE5DpjjAe4GXgB2Ao8bIzZHMl2hlKsHXOsHS/E5jFHo1j7HWPteCH2jjnWjjdaxeLvGGvHHGvHC7F5zAMRY0yk26CUUkoppZRSagQmTY+cUkoppZRSSimLBnJKKaWUUkopNcloIKeUUkoppZRSk4wGckoppZRSSik1yWggp5RSSimllFKTjAZySimllFJKKTXJaCCnlFJKKaWUUpOMBnJKKaWUUkopNcloIKeUUkoppZRSk8z/Bw7AAf3k0AahAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1080x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "all_metrics = state_history.unpack_stored_metrics()\n",
    "\n",
    "fig, (ax0, ax1, ax2) = plt.subplots(nrows=1, ncols=3, sharex=True, figsize=(15, 3))\n",
    "\n",
    "ax0.plot(all_metrics[\"iters\"], np.stack(all_metrics[\"eq_multipliers\"]))\n",
    "ax0.set_title(\"Multipliers\")\n",
    "\n",
    "ax1.plot(all_metrics[\"iters\"], np.stack(all_metrics[\"eq_defect\"]), alpha=0.6)\n",
    "ax1.axhline(0.0, c=\"gray\", alpha=0.2)\n",
    "ax1.set_title(\"Defects\")\n",
    "\n",
    "ax2.plot(all_metrics[\"iters\"], all_metrics[\"loss\"])\n",
    "ax2.set_title(\"Objective\")\n",
    "\n",
    "[_.semilogx() for _ in (ax0, ax1, ax2)]\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "5b5416066f703e73dc2bf3c4e66d036632ca71ae166bbdd9bd33a210b8a690df"
  },
  "kernelspec": {
   "display_name": "Python 3.8.12 ('cooper')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
